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A MATHEMATICAL ANALYSIS OF THE STABILITY 
OF A LUNAR FLYING VEHICLE WITH MANUAL CONTROLS 

By David Hall 


SUMMARY 


Equations of motion for a four-degree-of-freedom mathematical 
model of the lunar flying vehicle are developed and used to study the 
stability of the autonomous vehicle. The human pilot describing funa 
tion is used to represent the control system and derive the stability 
criteria for the closed-loop attitude control system. 


INTRODUCTION 


The system under investigation is a free-f light manually controlled 
vehicle operating in the lunar environment. A mathematical model of the 
vehicle dynamics and the pilot control function is used to study the 
system stability. The system is represented by the block diagram in 
figure 1. 

The vehicle consists of a platform which supports the pilot, pro- 
pellant tanks, engine, landing gear, and equipment. Representative 
configurations of thrust vector-controlled and kinesthetic-controlled 
vehicles will be discussed. The coordinates which describe the motion 
of the bodies are shown in figure 2. 

The vehicle dynamics are based on a four-degree-of-freedom mathe- 
matical model of plane motion, which considers two rigid bodies con- 
nected at a point which will be referred to as the gimbal point. Two 
coordinates X and Z locate the center of mass of one body in the 
plane. One coordinate 0 determines the angular position of the body 
with respect to vertical, and the fourth coordinate 0 determines the 
angular position of the second body with respect to vertical. 
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2 


For the purposes of this study, several assumptions were made which 
simplify the analysis but retain the fidelity of the model. 

1. Flight time is small (approximately 4 minutes). 

2. Mass, center of mass, and inertia of each body are constant. 

3. Pilot describing function is time invariant and can be repre- 
sented by the Laplace transform. 

4. Pilot control force is applied normal to the thrust vector 
which is along the coordinate axis of body 2. 

5. Pilot control force reaction is applied at the center of mass 

o, ** 1. ^ W, «i.. 

The equations of motion are developed using Lagrange’s equation 
and solved for the angular accelerations. The Jacobian of the solution 
allows the equations to be written in state variable matrix notation. 
Routh's criteria are used to determine the stability requirements for 
the vehicle characteristic equation. 

Vehicle control is represented by the experimentally determined 
human pilot describing function. The closed-loop system characteristic 
equation is derived and Hurwitz criteria used to study the effect on 
stability of varying the vehicle and control function parameters. 

The mathematical model has been programed for the digital computer 
and particular solutions have been generated to verify the analytical 
results. Although this model is not an analog of the manually con- 
trolled vehicle with variable mass characteristics , it will hopefully 
provide some insight into the stability and control of vehicles in this 
class. The philosophy is that an approximate model will do in the 
beginning and that we will enhance our theoretical understanding by 
attempting engineering applications and learning from them. 

Work related to this analysis is presented in the Bibliography. 


SYMBOLS 



matrix element 



coefficient of 


n 


s 


in characteristic equation 


damping coefficient 
dissipation function 
system error 

base of natural logarithms 
control force 

transfer function 
acceleration of gravity 
Hurwitz matrix 

principal minor of Hurwitz matrix 

moment of inertia 
Jacobian matrix 
pilot gain 

spring coefficient 

directed length from body 1 center of mass to gimbal 
directed length from gimbal to body 2 center of mass 
directed length from gimbal to application of control force 
mass 

generalized force 
generalized coordinate 
Laplace transform operator 
engine thrust 


kinetic energy function 
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general lag time constant 
general leadtime constant 

lag time constant approximation of the neuromuscular system 

»* 

t time 

V potential energy function 

X coordinate of body 1 center of mass 

x^ body axis 

Y controlled element transfer function 
c 

\ Y pilot describing function 

Z coordinate of body 1 center of mass 

body axis 

$ coordinate of body 2 angular position 

A determinant of matrix 

\ <5 gimbal angle 

1 damping ratio of the second-order component of the 

\ neuromuscular system 

0 coordinate of body 1 angular position 

! x pilot res^ction time lag 

aijy undamped natural frequency of the second-order component of 

I the neuromuscular system 

■< 

Subscripts; 

;j 1 body 1 

2 body 2 

f 

■V 
s. ; 
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b body 

c control 

i Input 

N neuromuscular system 

n integer 

p pilot 


Operators ; 

j / \ 

- derivative of ( ) with respect to 


a 


partial derivative of ( ) with respect to q 

( * ) derivative of ( } with respect to time 
(*' second derivative of ( ) with respect to time 

^ approximately equal to 


ANALYSIS 


Vehicle Equations of Motion 

The analysis of vehicle dynamics is limited to plane motion to al- 
low the effect of vehicle design on stability to be emphasized without 
unnecessary complexity. This results in a system with four degrees of 
freedom and consisting of two rigid bodies connected at a point. 

The vehicle equations of motion can be derived using Lagrange 1 s 
equation 



(X) 


0 

The functions are 


T = ~ M 1 (x 2 + Z 2 ) ' 1 ’ " 2 ' 1 


,[x 


p *2 P'P 

+ | i^e" + | mJ % + z tf + 


+ L 2 2 0 2 + 21^6 (X oos 0-2 sin 0) 

+ 2L 2 e(X cos g - Z sin g) + 2L L £ eg cos(9 - g) 


1 •? 

+ 2 V 


d = |c(e - g) 2 

V = M^gZ + Mgg/z + L x cos 9 + Lg cos gj + i- k(e - g) 2 

Let q.^ = 0 

*2 = 3 
q. 3 = X 

*4 = Z 

Then Q = -TI. sin(0 - 0) 

Q = -F L 
^2 c c 

Q = T sin $ + F cos 6 
3 c 


= T cos 3 - F c sin 9 


Differentiating, substituting into tbe Lagrange equation, and rearranging t^he equations into 
matrix form give the equations of motion 
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The solution of equation (2) will he of the form 

0 = f^e.e.g.e.Fj 
3 = f 2 (e,e,M,F c ) 

X = f 3 ^6,0,3,S,F o ) 

z = f 4 (e,0,e,6,F c ) 

Since the accelerations are independent of the linear position and 
velocity, the stability of the rotational motion may be investigated 
separately. Consequently, equation (2) was solved for the angular ac- 
celerations by premultiplying by the inverse of the coefficient matrix. 
Letting A represent the determinant of the coefficient matrix, the 
angular accelerations are 


(e - 3) 


-»0| 

- c ("i 4 “ 2 )j(“i * "s) 1 ; * “i " 2 L 2 [ L 2 4 \ “■<» - - s > 

- M 1 2 M 9 2 L 1 li L g 2 cos (0 - 2 )sin(e - g) 0 2 

- mfo + Ws 2 + 

- Vhihh + W 2 2 + Va) 3 ^ 6 - 

- M 2 F c L lj( M l + ^)[ J 2 - M 1 L 2 L c C0S < 6 - 6 >] ' 

+ M 1 M 2 L 2 2 sin 2 (0 - 3 )| ( 3 ) 
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AS - k(M x + + Mg)l^ + M^g^^ + Lg cos(9 - S)^j| ( 6 - g) 

+ °( M 1 + M s) J( M 1 + M 2 ) X 1 + Wx[ L l + L 2 °° S < 9 - 6) ]| ( ® " &) 
+ Wl^l 1 ! + M 1 M 2 L 1 2 + ^l) 3 ^ 9 ~ ^ 


+ “iVhV 00 ^ 9 - e)sin(6 - g)g ! 
+ M 1 2 MgTL 2 Lg cos(8 - 3)sin(0 - g) 


- ' M 1 + M 2 Pc Vl 


r 


L + L_ cos(0 

C u 


-*o 


+ M 1 L C J 1 + M 2 L 1 


(4) 


where 


A = (Ml + < 


Wl + M 1 L 1 2 ) + h*l(h * M 2 L 2 2 , 


+ M 1 2 Mg 2 L 1 2 Lg 2 sin 2 ( 0 - g) 


(5) 


Vehicle Stability 

Equilibrium is defined as that position of the variables for which 
the derivatives of the variables are zero. From equations (3), (4), and 
(5) 9 the equilibrium position for the rotational motion of the autono- 
mous vehicle is (0 - $) = 0. To study the vehicle stability in the 
neighborhood of equilibrium, the Jacobian is used to obtain a linear ap- 
proximation of the acceleration equations. 
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The Jacobian is 


!fi 

3f x 


3f l 

8q -l 

3q 2 

3<1 3 

3<1 4 


3f 2 

3f 2 

ffg. 

3q i 

3q 2 

^3 

3 % 

3f 3 

3f 3 

3f 3 


8 *1 

3q 2 

3q 3 

3< l4 

3f 4 

3f 4 

3f 4 

!£i 

3q l 

3q 2 

3q 3 

3< ii t 


where f = 0 



The equations for rotational motion may then be written in state variable 
matrix form 


• 

0 


0 1 C 0 
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0 

8 

= 

a 21 a 22 a 23 a 24 


• 

e 

+ 

a 25 

• 

6 


0 0 0 1 


0 


0 

• • 

0 


a 4l a 42 a 43 a 44 


• 

0 




( 6 ) 
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where 


a 2X 

a 22 = 
a 4l = 

a 42 = 


- a 23 = "M k ( M l + M 2)R M 1 +M 2 ) X 


+ M 1 TL 1 ( M 1 I 2 + M 1 M 2 L 2 2 + M 2 I 2 )} 


'2 + M 1 M 2 L 2( L 2 + L l ) 


" a 24 " - I C ( M l + M 2 )[( M X + M 2) X 2 + M X M 2 L 2 ( L 2 + L l)_ 


-a 


43 = I { k ( M l + M 2 )[( M X + M 2 ) I X + M X M 2 L x( L X + L 2 )] 


+ M xV L X ?L 2 


} 


a 44 = I c ( M x + m 2 )[( m x + m 2 ) i x + Wx( L x + hj 


( M x + m 2 ) 


Wx + M X L X 2 ) + h*l{h + M 2 L 2 2 


J 


The adequacy of the linear approximation of the equations of motion 
is illustrated by figures 3 and 4 which were obtained by numerical inte- 
gration of equations (2) and (6), respectively. 

The third column of [A] is the negative of the first column; conse- 
quently, the matrix is singular and the system has infinitely many 
equilibrium states. This agrees with the statement that the vehicle is 
in equilibrium for 0 = £. 


The vehicle characteristic equation is obtained by 

det£s[l] - [A]] = 0 


where [I] is the identity matrix and [A] is the coefficient matrix. 
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The characteristic equation is 


s + ( -a 


(- 


'22 " a Mt) s3 + ( a 22 a 44 " a 21 - a 2U a lt2 “ a 4s) ^ 
+ ( a 21 a 44 " a 24 a 4l + a )*3 a 22 “ a 23 a 42) S 
+ (%3 a 21 ” a 23 a 4l) = 0 


As a result of the relationships between the a.. given above, the last 

**• J 

two coefficients are zero, and part of the coefficient of s*~ cancels. 
The characteristic equation reduces to 


+ ( -a 22 - a 44) s3 + (" a 2 l - a 4 3 ) s2 = 0 


The stability of the vehicle may be investigated by Routh’s criteria. 
The Routhian array for a fourth-order system is 


b b 

, 0 2 
2 " * b 


b_ - 


b i b 4 


b t), 

*t_ o 3 

b 2 - ^ 7 “ 


b 2 b 4 
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w.iere the b ± refer to the coefficient of s x in the characteristic 

equation. Routh's criteria state that the number of roots of the char- 
acteristic equation with positive real parts is equal to the number of 
changes of sign of the terms in the first column of the Routhian array. 

For the vehicle characteristic equation 


b 

o 


= 1 


b 


3 



0 


and the Routhian array is 


0 


1 b. 


b 2 


The zero in the s row indicates that there are roots that are the 
negatives of each other. From the auxiliary equation 



= 0 


both roots are zero. There will be no changes of sign in the first col- 
umn and, therefore, no roots with positive real parts if 


\ > 0 


In terms of the vehicle parameters 


b l = " T U fr x M T 

I 2 M 2 , V I 1 + M 1 L 


+ m i m 2 ( l i + l 2 ) 2 1 

X 2 ) + h U l( X 2 + M 2 L 2 2 ) 


( 8 ) 


r" «■ 

k ( M l + M 2 )L( M 1 + ^(h + h) + M 1 M 2 ( L 1 + L 2 ) 2 _ 


b 2 " 


' a 21 ' a 43 


+ 

M n TL_ 

(m 

+ + M_ M-L_ / 

'l. + L 0 1 



1 1 _ 

A 1 

2/2 12 2' 

, 1 2) 


/V + M 

LM n (I 


M,L 2 J + I_M_ (I„ 

+ m 0 l 2> ) 


\ 1 2/ 

2 2 \ 

1 

11/ 11(2 

2 2 ) 



(9) 

Since all the terms in equation (8) are positive, b^ is always greater 

than zero. The denominator of equation (9) is .positive; therefore, the 
sign of b^ will depend on the sign and magnitude of the terms in the 

numerator. 

For example, with the gimbal point below the center of mass, 

is negative (fig. 2 and definition of symbols),, Assuming L is posi- 

et 

tive, the following equation must be satisfied for nondivergent oscil- 
lations . 
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In the kinesthetic-controlled vehicle concept, it is assumed that the 
pilot acts as a rigid body pivoted at his ankles. Therefore, the mass, 
inertia, and characteristic length of body 1 are fixed, and the design 
of the platform is constrained to be compatible with the stability cri- 
teria. 

The existence of a double root at s = 0 indicates that the origin 
is unstable. This can be understood by examining the block diagram of 
the vehicle equations of motion (fig. 5). The initial conditions 


9(0) = 6(0) + 0 


will produce a linear and equal increase in 6 and 3. 
For zero input, the solution of equation (6) is 


[d(s)] = [ S [I] - [A]]" 1 [d(0)] 


( 11 ) 


where 



is the Laplace transform of the fundamental matrix. 


For the vehicle system, the matrix is 


16 




IT 


The matrix shows that the solution will diverge for arbitrary ini- 
tial conditions. However, if the damping coefficient c is zero 


a 22 ~ a 24 = a l*2 



and the solution is stable for arbitrary displacements and zero veloc- 
ities. Therefore, misalinement between the thrust vector and the vehicle 
center of mass at lift-off is not critical since the oscillations do not 
diverge. However, any disturbance in flight will result in divergent 
oscillations unless corrected by the pilot. An example is shown in fig* 
ures 6 and 7 • 


The form of the response is determined by the denominator of equa- 
tion (12), which is the characteristic equation. Therefore, the natural 
frequency of oscillation is dependent on the location of the roots of 
the quadratic factor. Since the location of the gi.mbal point will ef- 
fect the value of (-a^ - a^), as shown in equations (9) and (10), it 

will effect the location of the roots. Obviously, the magnitude and 
frequency of the response is a function of gimbal point location. 


Pilot Describing Function 


The control system can be characterized as a manual process for 
minimizing visually perceived errors by exercising continuous control 
so as to match visually presented input and output signals. The pilot 
closes the control loop by determining the error between a reference 
attitude and the present attitude; this process is called compensatory 
tracking. 

The output of a nonlinear element (e.g., a human pilot) always con- 
tains higher harmonics, the effect of which depends on the characteris- 
tics of the other components. A quasi-linear equivalent to the nonlinear 
element, for a specific input, can be represented by a describing func- 
tion (the equivalent linear element) and a remnant (the difference be- 
tween the output of the nonlinear system and the equivalent linear 
element). Since components with relatively large inertia and long time 
constants filter higher harmonics considerably, a reasonable description 
for the pilot response characteristics is a linearized model of an 
equivalent transfer function which is linearly correlated with the sys- 
tem input. 
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To develop some criteria for stability of the controlled vehicle, 
it is necessary to have an expression for the control system. The pilot 
describing function has been used to express the functional relationship 
between the input and output of a single-axis manual control system 
(McRuer , et al . ) . 

For this analysis, the input is the pitch attitude error 



( 13 ) 


and the output is the pilot limb position (i.e., the angle between 
body 1 and body 2; the gimbal angle) 


5 = 6-3 


( 14 ) 


Thus, angular control of the vehicle in the pitch plane is represented 
by 



Since the parameters of the describing function are adjustable by the 
pilot, the function is only valid in the frequency domain and only ex- 
ists under essentially stationary conditions. 


Although the describing function represents a functional relation- 
ship, there appears to be some structural correlation. For example, the 

«7'C2 

time delay e k is due to excitation of the retina, nerve conduction, 
and data processing in the central nervous system. The neuromuscular 
system characteristic 

\ 



2 ? n s 


03 


N 
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i 3 compatible with physiological data for the muscle spindle and mu:,.cle 
tension, and the behavioral characteristics of the muscle /manipulator 
combination, The manipulator or control stick dynamics are not explic- 
itly present in the term since virtually all experiments have been con- 
ducted with manipulators of negligible inertia, and the Information which 
exists is insufficient to account for the effect analytically. However, 
the manipulator dynamics can be represented in terms of the gimbal angle 
and control force, using equations (6) and (l4), as 

& + (~ a 22 " + (” a 21 " a l * 3 ) ^ = ( S 25 _ % 5 ) F c 


and in Laplace transform notation as 



The form of the manipulator dynamics is the same as the second-order 
component of the neuromuscular system. Since existing data support a 
first-order approximation to the neuromuscular system, the system is 
usually represented by T^s + 1 where 




25 


N 


00 . 


N 


The factor 


T t s + 1 
Jj 

TjS + 1 

in a compensation or equalizing characteristic which, together with the 
gain K , is adjustable by the pilot to obtain system stability. 
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To simplify the subsequent analysis, the time delay in equation (15) 
can he replaced by a Pade approximation 





The pilot describing function then becomes 


Y 

P 


K P (v + l )(- i s + 1 ) 

e ( T i s + 1 )(v + 1 )(i s + 1 ) 


(17) 


Attitude Control System Stability 

From equations (6) and (l4), the linear equation for vehicle pitch 
attitude 0 can be written in terms of the gimbal angle 6 as 


0 = a 2] _fi + a 22 6 


( 18 ) 


Taking the Laplace transform and rearranging, the transfer function for 
the vehicle dynamics is 


Y 


c 


_0 

5 


a 22 S + a 21 


(19) 


The product of the transfer function for the control system (eq. (17)) 
and the transfer function for the vehicle dynamics yields the system 
open-loop transfer function 


G(s) = Y X 
P c 



(t l s + l)(- f s + l)(a 22 

s 2 (v + i )(v + l )(l s 


S + a 2l) 


( 20 ) 


and the closed-loop transfer function is 
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where 



A polynomial whose zeros all have negative real parts is called a 
Hurwitz polynomial. The necessary and sufficient condition for equa- • 
tion (22) to be a Hurwitz polynomial is that the principal minors of 
H all be positive. 
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The principal minors of H are 




From equation (23) it can be seen that 


b > 0 
o 

\ > 0 


since all the parameters are positive. 
Therefore 



> 0 


i 


remaining principal minors, vhich must be positive, are 


2k 
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The equations for the principal minors indicate that Hg, H^, and 

Hj. could he positive or negative depending on the sign and magnitude of 
an< ^ a 22* parameters T^ s T^, and must he adjusted by 

the pilot to insure that the minors are positive, thus rendering the 
system stable. For a particular set of vehicle parameters and values of 
T^ and t, the inequalities for H n may he solved to obtain the per- 
missible range of values of T^, T^ , and K . 

The impulse response of the attitude control system closed-loop 
transfer function is shown in figure 8 for a set of parameters which 
satisfy Hurwitz criteria and in figure 9 for a set which do not. 


CONCLUSIONS 


The autonomous vehicle is unstable except for the special case of 
zero damping and initial conditions on displacement only. The magnitude 
and frequency of the response differ for the configurations with the 
gimbal point above or below the center of mass since the vehicle param- 
eters which define the configuration determine the location of the roots 
of the characteristic equation. 

Kinesthetic control allows less freedom in design because the phys- 
ical parameters of the pilot (i.e., mass, inertia, length) cannot be ad- 
justed to meet the stability criteria. 

The stability criteria for the vehicle and for the closed-loop sys- 
tem should be useful, in defining the permissible range of vehicle param- 
eters for stable operation. Extension of this procedure to multiaxis 
control will be complicated by vehicle axis coupling through products of 
inertia an$ by the requirement to generate different equilizations in 
different axes simultaneously. 
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Figure 1.- Single-axis attitude control system block diagram. 




Figure 2.- Coordinate system 
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Figure 3.- Numerical solution of nonlinear equations of motion. 
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Figure 6.- Response of thrust vector controlled vehicle (autonomous, 

gimbal above center of mass). 
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Figure 8.- Impulse response for attitude control system (stable) 
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Figure 9«- Impulse response for attitude control system (unstable) 









